clear, clc;

% 城市地图
C = zeros(,2)
% M为城市数量，N为矩阵维度
[M,N] = size(C)
% 距离矩阵
distance = zeros(M,M);
% 求距离
for m = 1:M
    for n = 1:M
        distance(m,n) = sqrt(sum((C(m,:) - C(n,:)^2)));
    end
end

m = 50;    % 蚂蚁个数
G = 150; % 迭代次数
% 信息素参数 
alpha = 1; % 信息素重要程度（一般1~4）
beta = 5;  %启发式因子的重要程度（一般3~5）
rho = 0.25; % 信息素蒸发系数
Q = 100; % 信息素增加系数
Eta = 1./distance; % 启发式因子
Tau = ones(M,M);  % 路径的信息素矩阵
% 路径信息
Tabu = zeros(m,M); % 禁忌表（记录每只蚂蚁走过的路程）
gen = 1;
R_best = zeros(G,M); % 各代的最佳路线
L_best = inf.*ones(G,1); % 每一代的最佳路径的长度（初始设为无穷大）

% 开始迭代计算
while gen < G
    random_pos = [];
    % m只蚂蚁放到M座城市
    for i = 1:(ceil(m/M)) 
        % 每只蚂蚁放到随机的城市
        random_pos = [random_pos, randperm(M)];
    end
    Tabu(:,1) = random_pos(1,1:m)'; % 第一次迭代中，每只蚂蚁的禁忌表

    % 从第二个城市开始循环
    for i = 2:M
        % 每只蚂蚁的循环
        for j = 1:m 
            visited = Tabu(j, 1:(i-1)); % 在访问第i个城市时，第j个蚂蚁访问过的城市
            unvisited = zeros(1,(M+1-i)); % 待访问的城市
            visit_P = unvisited; % 蚂蚁j访问剩下的城市的概率
            count = 1;
            % 找出未访问的城市
            for k = 1:M
                % 判断该城市是否还没有访问过
                if isempty(find (visited == k)) 
                    unvisited(count) = k;
                    count = count + 1;
                end
            end
            if gen >= 2
                Tab(1,:) = R_best(gen-1,:);
            end
            % 记录m只蚂蚁迭代的最佳路线
            L = zeros(1,m);
            for i = 1:m
                R = Tabu(i,:);
                % 走一周回到起点
                L(i) = distance(R(M),R(1)); 
                for j = 1:(M-1)
                    L(i) = L(i) + distance(R(j), R(j+1));
                end
            end
            % 记录每次迭代中路径最短值
            L_best(gen) = min(L);
            pos = find(L==L_best(gen));
            R_best(gen,:) = Tabu(pos(1), :); % 最优路径

            % 更新信息素的值
            Delta_Tau = zeros(M,M);
            % m只蚂蚁遍历
            for i = 1:m
                % M个城市遍历
                for j = 1:(M-1)
                    Delta_Tau(Tabu(i,j),Tabu(i,j+1)) = Delta_Tau(Tabu(i,j),Tabu(i,j+1)) + Q/L(i); % m只蚂蚁的信息素累加 这里采用的是论文中ant-cycle模型
                end
                Delta_Tau(Tabu(i,M),Tabu(i,1)) = Delta_Tau(Tabu(i,M),Tabu(i,1)) + Q/L(i);
            end
            % 更新路径上的信息素含量
            Tau = (1-rho).*Tau+Delta_Tau;
            % 禁忌表清零
            Tabu = zeros(m,M);

            for i=1:(M-1)
                plot([C(R_best(gen,i),1),C(R_best(gen,i+1),1)],[C(R_best(gen,i),2),C(R_best(gen,i+1),2)],'bo-');
                hold on;
            end
            plot([C(R_best(gen,n),1),C(R_best(gen,1),1)],[C(R_best(gen,n),2),C(R_best(gen,1),2)],'ro-');
            title(['最短路径：',num2str(L_best(gen))]);
            hold off;
            pause(0.05);
            gen = gen+1;
        end
        figure(2);
        plot(L_best);
